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I. INTRODUCTION 


With the high accuracy of guidance information provided by the Microwave 
Landing System (MLS), the use of curved flight paths and steep glideslopes in 
the final approach in terminal areas has become possible. As the Microwave 
Landing System is less sensitive to weather conditions than conventional systems, 
automatic landing systems using the MLS can be used to reduce the congestion 
in terminal areas due to adverse weather conditions. To meet projected aviation 
requirements through optimum utilization of the MLS and "state of the art" 
guidance and control technology, NASA has introduced the Terminal Configured 
Vehicle (TCV) Program at the Langley Research Center, [l]. The goals of the 
TCV Program include the reduction of aircraft noise in airport communities, all- 
weather operation in the, terminal area, and the efficient utilization of airspace in 
congested terminal areas. 

The Microwave Landing System is a ground-based guidance system under 
development by the FA A, which provides position information to aircraft within 
its volumetric coverage [2]. It consists of a DME providing range information, 
an azimuth antenna generally collocated with the DME antenna providing the air- 
craft's azimuth angle relative to the runway centerline up to + 60° , and an 
elevation antenna located at the glidepath intercept point but offset to the side of 
the runway providing the aircraft's elevation angle up to 20°. The MLS thus has 
a volumetric coverage, and provides guidance information that can be used for 
steep approaches and curved flight paths in the terminal area. 

This work considers the development of a digital automatic control law for 
the TCV research aircraft (B-737) to perform a steep final approach in automatic 
landings. The use of a steep glideslope in the final approach reduces the intensity 
of the noise perceived on the ground, provides a method to avoid the vortex gener- 
ated by a large aircraft flying a shallower glideslope, and provides more flexibility 
in sequencing and spacing of aircraft at different altitudes, thus allowing more 
efficient use of airspace. A 3D, constant airspeed control law which engages at 
localizer capture and brings the aircraft to touchdown was developed for the lateral 
and longitudinal dynamics. Localizer and glideslope capture can be engaged simul- 
taneously, resulting in a curved capture path horizontally as well as vertically. 



A decrab law with a maximum roll of 5® is used to control the effects of cross- 
winds in landings. The flare law was developed to follow a fixed flare path in 
order to reduce touchdown dispersion under varying wind conditions in steep 
approaches. The control system uses MLS position data and body-mounted 
accelerometer data in place of inertial platform accelerations, which are filtered 
along with attitude and body rate information. The filter provides less noisy 
estimates of the data obtained by the instruments as well as estimates of non- 
measured parameters such as the sideslip angle and wind velocities. These 
filtered outputs are then used to control the aircraft during the final approach. 

The control law considers the localizer capture, glideslope capture, steep glide- 
slope tracking, and flare to touchdovm phases of flight and includes a decrab law 
for automatic landings where a crosswind is present. The control law was devel- 
oped using modem digital optimal control techniques which are well-suited to the 
discrete character of the Microwave Landing System. 

In Section II, the aircraft’s equations of motion and the wind model are des- 
cribed. A perturbation model of the aircraft’s lateral dynamics for a steep 
glideslope of 6^ is obtained and is expressed in state variable form. Using the 
Dryden spectrum, a dynamical model for the simulation of wind gusts affecting 
the aircraft’s lateral dynamics is developed. Then steady winds and shear are 
included into the wind model. 

In Section III, a mathematical model for the noises in the various on-board 
sensors and a model for the noisy MLS data are developed for simulation purposes. 
The development of the filter equations is then described. 

In Section IV, the development of a digital automatic control law for the 
lateral dynamics of the aircraft and a constant airspeed, 3D control law for the 
longitudinal dynamics is described. The control uses the filtered estimates of 
measured variables as well as wind velocity estimates to control the aircraft’s 
deviations from runway centerline and the glideslope during the final approach. 

In Section V, the results obtained from a digital simulation of the aircraft 
dynamics, wind conditions, sensor noises using the control law and filter devel- 
oped are described. 
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IL MODEL OF THE AIRCRAFT DYNAMICS 
AND WIND CONDITIONS 

In this section, mathematical models for the lateral and longitudinal 
(including vertical) dynamics of the aircraft and the wind components which have 
an effect on the aircraft’s dynamics will be developed. The models developed 
will be used in the digital simulation of the motion of the aircraft under different 
wind conditions and in the development of the filtering algorithm and the automatic 
control law. The equations thus obtained will then be expressed in state variable 
form which is more suitable for the application of modem optimal control tech- 
niques. The control law developed will then be evaluated using a non-linear model 
of the aircraft d3mamics. 

A. MODEL OF THE LATERAL DYNAMICS 

To simulate the motion of the aircraft on a digital computer, a mathematical 
model is necessary. The general equations of motion for aircraft are non-linear 
differential equations where the longitudinal and lateral motion variables are 
coupled [3] , [4] , [5] . However, under a steady flight condition assumption, 
these equations can be simplified significantly with little change in their validity 
[3, pp. 154-165], [4, p. 2.30] . In the development presented here, the phases 
of the final approach considered are the localizer capture, glideslope capture, 
glideslope tracking and flare. With the exception of localizer capture, in the 
phases of flight considered the roll, yaw and pitch angles have quite small values, 
where the yaw angle is measured relative to runway centerline; thus, with a con- 
stant airspeed, a steady flight condition would be present. During localizer cap- 
ture, the yaw angle may have large values. Similarly, the roll angle may also 
have large values during capture; however, as the time period over which large 
values in the roll angle would be sustained is relatively short, this would not 
ihvalidate the simplified equations of motion. On the other hand, the yaw angle 
does not affect the values of the stability derivatives in aiygiven flight condition. 
This is due to the fact that roughly the yaw angle describes the angular deviation 
of the body x-axis from the runway centerline, whereas the stability derivatives 
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depend on the values of aerodynamic variables such as the sideslip angle, angle 
of attack, etc, , which describe the motion of the aircraft relative to the atmos- 
phere surrounding it. The yaw angle is used, however, in the computations of 
the aircraft’s position and velocity relative to an earth-fixed coordinate frame; 
i. e. , in transforming body-axis variables into inertial variables. Thus, for 
large values of yaw, the small angle approximation in the transformation from 
say stability axes to earth-fixed axes does not hold and these transformations 
have to be dealt with separately. 

When the steady-state body attitude rates, the steady-state velocity com- 
ponent along the body y-axis and the steady-state bank angle have a value of 
zero, the general equations of motion can be linearized to obtain the perturbation 
equations, where the longitudinal and lateral variables are no longer coupled 


, see [ 2 ] , p. 2.33). 



m ( V + Uo r - Wo p ) 

= mg cosyo 9 + fAy + % 

(1) 

if 

11 

+ -t/j- , 

(2) 

Izz ^xz P ~ ”A 

+ n^p 

(3 ) 


Equations ( 1 ) - ( 3 ) are the perturbation equations which describe the lateral 
dynamics of the aircraft, where 


Uo 

Wo 

V 

P 

r 

Vo 

9 

^Ay 

% 


steady-state inertial speed along the stability x-axis, i. e. along Xg , 

steady-state inertial speed along the stability z-axis, i. e. along Zg , 

perturbation in the inertial speed along the stability y-axis, i. e. 
along yg , 

perturbation in the roll rate about Xg , 
perturbation in the yaw rate about Zg , 
steady-state flight path angle, 
perturbation in the roll angle, 

perturbation in the net aerodynamic force along yg , 
perturbation in the thrust force along yg , 
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perturbation in the rolling moment due to aerodynamic forces, 
perturbation in the yawing moment due to aerodynamic forces, 
perturbation in the rolling moment due to thrust, 
nrp - perturbation in the yawing moment due to thrust, 

and , I 22 and are the moments of inertia of the aircraft in the stability axes. 

Now, the terms f^y , f'py > can be expressed in terms of 

the stability derivatives of the aircraft evaluated at the steady state values of 
the aerodynamic variables and the control surface settings, in linear form [ 3 ] , 

C 4 ] , C 5 ], Substituting these expressions into equations ( 1 ) - ( 3 ) and 
rearranging terms we obtain linear differential equations in the sideslip angle, 
the roll rate and yaw rate. Writing the derivatives of the roll and yaw angles 
in terms of roll and yaw rates, we obtain the following set of differential 
equations. 


9 = sec 0(3 (cos Vo p + cosy 0 r) (4) 

• 

If = sec6o (sin a© P + r) ( 5 ) 

P = 9 + S33 3 + 834 p + 33 5 r + bai 6 A + b 32 6 R 

+ bss 6 sp + + hsa p^ + II33 r^ ( ^ ) 

p = 843 3 + 844 p + 845 r + b4i 6 a + b42 6 r + b43 5 sp 

+ h4rPw*^^Pw + h43iNv (7) 

r = 853 3 + 854 P + ^55 ^ + bs3 6 r + 1%3 6sp 

+ h 5 ^Pvv'’"^Pw'^^^w (8) 


where 3 , P and r are the sideslip angle, roll rate and yaw rate relative to the 
inertial frame of reference, respectively, 3^ t Pw ^w nre the sideslip angle, 
the roll rate and yaw rate due to wind velocities only, 6 a, 6r and 6sp are the 
perturbations of the ailerons, rudder and spoilers, respectively. The coefficients 
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ajj in the above equations depend on the aircraft stability derivatives and are 
given in the Appendix. Thus, a set of linear differential equations describing 
the lateral velocities and attitude of the aircraft can be obtained. 

The position of the aircraft relative to runway centerline is expressed by 
the perpendicular distance of the aircraft center of mass to the runway center- 
line, This distance normalized by the aircraft’s steady state speed will be 
used as a state variable in addition to the equations already obtained. Thus, 
let Les be the matrix representing the transformation from the aircraft stabil- 
ity axes to the earth-fixed axes; so that Lgs (i, j ) is the element in the i*^ row 
and jth column of the matrix. Then, the lateral distance y (in feet) of the air- 
craft can be expressed by 

y = Uo[ Les ( 2, 1 ) ( 1 + u' ) + Les (2, 2 ) 3 + Les ( 2, 3 ) a] , ( 9 ) 

where u' is the normalized inertial perturbation in the speed along the Xg 
direction and a is the perturbation in the inertial angle of attack, and 

Les ( 2, 1 ) = cos o/q cos 6 sin i|r + cos cpsin sin 6 sin ijr 

- sin oto sin qjcos ijr ( 10 ) 

Les ( 2, 2 ) “ sin cpsin 6 sin + cos cp cos 'If (11) 

Les ( 2, 3 ) = - sin at, cos 6 sin + cos Oq cos cpsin 6 sin ijr 

- cos at) sin cp cos t . ( 12 ) 

Rearranging the terms in equation ( 9 ) 

yi = P + cos( ao - 6o) 'I' - sinao cp + Tly , ( 13 ) 

where T|y and y' are given by 

lly = ( Lgg (2, 2)-l)3 + Les (2, 1 ) ( 1 + u' ) - cos(ao - 6o) tlf 

+ Les (2, 3 ) a + sinoo cp , ( 14 ) 

y' = y/Uo . ( ) 
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Note that equation ( 13 ) contains no approximation when a is interpreted 
as the normalized inertial velocity component in the Zg direction, but is simply 
a rearranged form of { 9 ) with the non-linear terms grouped into a single term. 

The controls which affect the lateral motion of the aircraft are the aileron, 
rudder and spoiler surface settings as can be seen from equations ( 4 ) - ( 8 ) , 
where the spoiler action is used only to aid the effect of the ailerons during 
turns. Thus, the spoiler setting is programmed according to the aileron setting 
and is modelled here as 

6sp = Cgpg , Cgp3 = 1. 73 . { ) 

To further account for the relatively slow motion of the rudder, a rate 
command is preferred to a position command; hence, the rudder position is con- 
sidered a state variable which is obtained by integrating the rudder rate command. 

6r = ua , ( 17 ) 

where ug is considered to be the rudder rate control. 

If the relation between the spoiler and aileron given in ( 16 ) is substituted 
into ( 6 ), ( 7 ) and ( 8 ) , then the spoiler terms are eliminated from the equations. 
Now, forming a state vector x such that 

- ( cp ^ p p r y» 6 R ) ( 18 ) 

and a control vector u such that 

= ( 6 A 6 R ) , ( 19 ) 

equations ( 4 ) - ( 8 ) , ( 13 ) , ( 16 ) and ( 17 ) can be combined into a state 
variable model of the lateral motion of the aircraft of the form, 

X = Ax + Bu + Dw + 1] , ( 20 ) 

-T 

where Pw T|=(OOOOOTlyO). 
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Thus, equation ( 20 ) describes the lateral motion of the aircraft when the 
aircraft is controlled by the surface settings of the ailerons and rudder in the 
state variable format. The effect of winds on the motion is also included 
through the term w. 

B. MODEL OF LONGITUDINAL DYNAMICS 


The longitudinal equations, which describe the motion of the aircraft in the 
forward and vertical directions, can be obtained using a similar treatment as the 
lateral equations of motion. As noted in the last section, for the phases of flight 
considered in this work, the equations of motion can be described as perturbations 
about a steady flight condition where the longitudinal and lateral variables are 
uncoupled. Further, the pitch angle and the angle of attack are usually small 
enough to allow the use of small angle approximations. Thus, the longitudinal 
equations of motion can be described by linear differential equations in the per- 
turbations from the steady flight condition as shown below C 3 ] . 

m u = mg cosYo ©p + Qo S [ ( - Cd^, + 2 CDq + Ctxu' + 2 Ctxq ) 

u' + (Clo - a - CDfe Se-CD6g + Cxx^t ] (21) 

m ( w - Uo q ) = -mg sinvo 6p+ Qo S (Clui + 2 Clq ) u’ 

- ( + Cdo ) « - Ci^ a - CLq fl ] ( 22 ) 

lyy q = Qo S ( Cmu' 2 Cj^o ) H! + ( * CmT^) 

+ k + C]viq a_ + Cm^ 6e + Cmss ^M&r ^ 


ii! = 

u* + u’ 

w ’ 

II 

It 

( 24 ) 

or = 

^W ’ 

w 

""“Uo ’ 

( 25 ) 

a = 

q + q\y » 


( 26 ) 

®p = 

0 
CD 

1 

CD 


( 27 ) 
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u and a are the perturbations in the inertial speed along the x stability axis and the 

inertial angle of attack, respectively, q is the pitch rate, and are the 

components due to winds, qo is the steady state value of the dynamic pressure, 

S is the effective wing area, 6e, 6s and 6T are the perturbations in the elevator, 
stabilizer and thrust from their steady values, respectively. A more detailed 
derivation of the longitudinal equations of motion used here is given in [ 10 ] , 

[11]. 

The position of the aircraft relative to the earth-fixed coordinate frame can 
be obtained by integrating the inertial velocity components along each coordinate; 
thus, 

♦ 

X' = -j^= Les( 1. 1) (1+u') +Les(1. 2)P+Les(1. 3)0 , (28) 

Uo 

z’ = •^— = LeS ( 3, 1 ) ( 1 + u') + Les ( 3, 2 ) 3 + Les ( 3, 3 ) a . ( 29 ) 

Uo 

Note that oand 3 are inertial quantities and correspond to normalized velocity 
components in the stability axes; under no wind conditions these would be the same 
as the aerodynamic angle of attack and sideslip. We can rewrite the position 
equations in the following form. 

x' = - sinvo 6p + cosYo u' + sinYo o+ llx 

z* = - cosYo ®p “ sinYo u' + cosYo T|z 

where 

\ = Les ( 1 ) ( 1 + u* ) - cosYo u' + SinYo 6p + Les (1.2) 3 

+ ( Les ( 1, 3 ) - SinYo ) ot , ( 32 ) 

flz = Les ( 3. 1 ) ( 1 + u' ) + SinYo u' + cosYo 6p + Les (3,2) 3 

+ ( Les (3,3) - cosYo ) • ( 33 ) 

In this form, the differential equations are linear with a forcing function 
that contains the non-linear part which are second order terms with respect to 
the steady flight condition considered. With this approach we can use linear 
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(30) 

(31) 



theory in the development of the filter and control law without neglecting the 
non-linear terms completely. 

To account for the effects of the servo responses for the stabilizer and 
engine dynamics, simple linear models were used. 


6 T = - . 5 6 t + .298 6 th , 

(34) 

6 th = U 3 , 

(35) 

6 s = U 2 , 

(36) 


where 6 T is the thrust perturbation in units of one thousand pounds per unit 
of 6 T , 6 th is the throttle perturbation in degrees and 6 s is the stabilizer, 
pertuibation in radians. As the lags in the elevator action are small, the 
elevator time constant was neglected. The aircraft’s longitudinal equations 
of motion, the position equations and the actuator equations can be combined 
and after some manipulation can be expressed in state variable form. 


11 

CD 

O' q x’ z’ 

6 T 6 th 6 s) , 

II 

( 6 e 6 s 

6 th) , 


II 

( u'w “w 

<4w ) * 


11 

(0 0 0 

0 ^ Tiz 

0 0 0 ), 


+ Biui ^ \ , ( 37 ) 

where A^, B^, and D^are matrices of appropriate size corresponding to the 
coefficients in the original equations. Expressions for the elements of these 
matrices can be obtained in terms of the stability derivatives of the aircraft. 
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C. MODEL OF THE WIND CONDITIONS 


The lateral motion of the aircraft is described by the state variable model 
given in equation ( 20 ); this model describes the response of the aircraft when 
a control is applied or when the wind velocities such as gusts or steady winds 
are non-zero. The effects of the wind velocities are introduced through the 
vector w. The components of this vector are ®r the wind velocitiy along the 
yg direction normalized by the airspeed of the aircraft, p^^ or the rotation of 
the air around the aircraft about the Xg axis, and r^ or the rotation of the air 
around the aircraft about the Zg axis, respectively. The roll rate p^ and yaw 
rate r^ components of the wind vector w consist only of the effects of wind gusts, 
thus having an average value of zero; i. e. , these components do not have a 
steady state effect but introduce turbulence effects into the equations. On the 
other hand, the or the normalized lateral wind velocity contains terms for 
both wind gusts and steady winds; thus, 


Pw “ 3g + 3s > 


(38 ) 


where 3g is the gust or turbulence term, and pg is the steady wind term. The 
gust terms are of a random nature and can be modelled using the well-known 
Dryden spectrum [3 ] , C 4 ] . This method consists of using spectral factor- 
ization methods to obtain a dynamical system which generates a random process 
having the specified power spectral density when driven by a white noise process 
[6 ], [7 ], [8 ]. 


The Dryden spectra describe the statistical behaviour of wind gust veloci- 
ties along the aircraft body coordinates by specifying their power spectral densi- 
ties in terms of the spatial frequency 0 [ 9 ] . The spectra for the gust components 
of interest are given below. 


^ 2 Lv 1 + 3 

tpp (Cl) - Uy — ■ ■■ 

ttv^ Ci + (Lvnf]® 


( 39 1 
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«pp 


2 


( 40 ) 


Pw 

Lw 


g (TT_^3 

4b 

1 +(i^) 

TT 


Vr ( ^ ' 


1 +( 


3b fl y 2 
n 


(pp ( Q) 


( 41 ) 


where b is the wing span, Ly and are the scales of turbulence, V is the 
airspeed, a is the variance of the lateral gust and is the variance of 

the vertical gust. The change from the spatial frequency 0 to the temporal 
frequency cu can be made by 


n = 


UJ 

v 


( 42 ) 


Substituting equation ( 42 ) into the expressions for the power spectral 
densities of the wind gusts we obtain the spectra in terms of the temporal 
frequency; then using spectral factorization techniques the following transfer 
functions can be obtained 


H 3 (s ) 


Hp ( s ) 



1 + ^ ^ ® 
V 


V ' 

% 

r 1 -]h 



Ll^ vJ 

1 



TT V 


( 43 ) 


( 44 ) 


( s ) 


1 + 


3b s 
n V 


H 3 ( s ) 


( 45 ) 


It should be noted that even though p^ is independent of and , the 
latter two are not independent of each other Thus, if a white noise process 
is input to the transfer function Hj. ( s ), the output would have the desired 
power spectral density, but may not have the desired cross-correlation with 
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pg . Hence, equation ( 45 ) must be interpreted as Pg being the input to 
the first term in the above equation in order that the proper cross-correlation 
be obtained. 

The wind gust terms can thus be simulated by passing white noise 
through the systems with the transfer functions given in equations ( 43 ), ( 44 ), 
and ( 45 ). The lateral wind, however, has a steady or average value which 
is not necessarily negligible. Thus, consider that a steady wind is present; 
in the earth-fixed coordinate system, the wind velocity has a component in 
the direction of runway centerline Wjj , and a component perpendicular to the 
runway centerline say Wg; it is assumed that there is no steady wind in the 
vertical direction although gusts may be present. Hence, if L 5 £ is the 
transformation matrix from earth-fixed to stability coordinates, then the 
steady component, Pg , of the normalized lateral wind velocity is given by 

Ps = LsE ( 2. 2 ) Ws + Lse ( 2, 1 ) Wx , ( 46 ) 

Lsj; ( 2, 1 ) = sin 0 sin tp cos i|r - cos <p sin i|; , ( 47 ) 

Lgg ( 2, 2 ) = sin 6 sin {p sin ^ + cos cos i|f . ( 48 ) 

To include wind shear into the model, the steady lateral wind velocity 
can be described as 

Wg = Wa + 0)3 , ( 49 ) 

W = u >4 , ( 50 ) 

6 

where 0)3 and 0)4 are gaussian white noise processes independent of each other 
and of Pg , P;y and . Now, the transfer functions for the gusts described 
in ( 43 ), ( 44 ), and ( 45 ) can be combined into a state variable model of 
fourth order. Adding ( 49 ) and ( 50 ) to this model we obtain a sixth order 
model of the form 

W = W + u) , ( 51 ) 
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where is Pg , W 3 is , W 4 is , u> is a 4-vector of independent 
gaussian white noise processes, and are given in the Appendix. 

Thus, the wind vector W can be generated using equation ( 51 ). To 
complete the development of wind conditions the vector w is needed for sub- 
stitution into ( 20 ); this vector can be obtained from W as follows. First 
note that the elements of w are in the stability axes; hence, Wa and W 4 must 
be transformed from body to stability axes; then the steady winds must also 
be expressed in this coordinate system. Thus, we obtain 


Pw 

= + W5 + , 


Pw 

= sin (Xq 

W3 + cos Ofo 

W4 , 


= cos ao 

W3 - sin Oq 

W4 , 

where § 

1 

= ( ^sE 

(2, 2) - 1 ) 

Ws + Lse (2 . 1) Wx 


These equations can be expressed in matrix notation as 
w = W + § , 


( 52 ) 
( 53 ) 
( 54 ) 

( 55 ) 


where 5 is a vector with its first element § as given above and its other 

1 

two elements zero. Thus W can be generated using ( 51 ), and w given in 
( 55 ) can be substituted into ( 20 ) to simulate the effects of a given wind 
condition on the motion of the aircraft; hence, 


X = Ax + Bu + DCw W + T1,71=D§+11 , (56) 

W = W + Bw «> . ( 57 ) 

The longitudinal wind model contains the components of steady wind velocities, 
turbulence and shear winds in the longitudinal axis. The turbulence model uses 
the Dryden spectra [ 4 ] for the various components varying with altitude. The 
turbulence model has three components: u’g in the direction » Q^g in the 
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directioiv and Qg which models the effect of turbulence on the pitch rate of the 


aircraft. 

These components are modelled using the following spectra. 


Su 

(G) = 

2<^u 

( 58 ) 

l + (LuG)^ 


(G) = 

1 + 3 (L^fl) 

( 59 ) 

“Va^ Ll + (LwG)^J^ 

Sq 

(G) = 

^ Va^ ^ 

, . , 4b ^ ’ 

(60) 


where b is the wing span, Lu and are the scales of turbulence, Va is the air- 
speed, and Qis the spatial frequency related to the temporal frequency ui by 

n = u)/Va . ( 61 ) 

The u'g component is independent of Og and qg, however, Og and qg are cor- 
related with their cross spectral density being 

Sqa («)) = S^(tu) . (62) 


The above spectra can be factored using spectral factorization methods to obtain 
a linear system driven by white noise which generates an output having the above 
spectral characteristics [6], [7], [s]. Thus, the following transfer functions 
are obtained to generate u’g , Og and qg . 


Gu (s) = 


Ga(s) = 


* ' 




Gq (s) 


, 4b 
1 + -=^ s 


TTVa 


(63 ) 


(64) 


(65) 
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where ofg is the input to the system Gq ( s ) to obtain qg with the specified spec- 
trum and cross-spectral density. Figures 1 and 2 show a block diagram of the 
system generating the turbulence components. 

The steady and shear wind in the longitudinal direction was modelled by 

^'s ~ ^'sh » u'sh “ t 

^S = ^4 • 

Thus, to simulate a specified shear profile for u'g , with appropriate 
initial conditions, e. g. , to obtain a linear profile u'g changing at a rate of u's^q , 
the initial condition for u'gjj is set to u'gjio and §^3 is set equal to zero; 
alternately, an impulse in will also achieve the same profile. 

The transfer functions obtained for the wind model can equivalently be 
expressed as differential equations in state variable form as shown in ( 68 ) ; 
the matrices B^^and C^^are given in the Appendix. 

Wi = Awi We, + It . wt = C^i V/t , (68) 

where w,ti® given by equation ( 37 ). 


( 66 ) 

(67) 
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nL DEVELOPMENT OF FILTER EQUATIONS 


This section describes the development of the filtering equations which are 
used with the longitudinal and lateral control laws during the final approach. The 
filters are used to reduce the noise in the various sensors and to obtain estimates 
of parameters which are not directly measured such as the wind velocities. As the 
Microwave Landing System ( MLS ) provides high accuracy position data to the air- 
craft at discrete points in time, a digital design was considered appropriate for 
the filtering function; this is also suitable for implementation on a digital computer 
such as the mini-computers on board the TCV aircraft. Furthermore, the perfor- 
mance of an analog ( or continuous ) filter implemented on a digital computer 
generally degrades as the sampling period or the integration step size increases , 
thus requiring a relatively high update rate, and computation time. The use of 
a digital design allows the choice of error feedback gains which take the sampling 
period into account as a constraint; thus reducing the degradation in performance. 
The form of filter used was a Kalman filter with constant gain: the steady state 
Kalman gains were used as a starting point to adjust the gains for good performance. 
A. MLS AND SENSOR MODELS 

The aircraft’s position is obtained using the MLS guidance system, which 
provides volumetric coverage in the terminal area. The aircraft receives range, 
azimuth and elevation information at discrete intervals from which it can obtain 
its position with high accuracy even under adverse weather conditions. The 
ground azimuth antenna is located at the runway centerline past the end of the 
runway; thus it provides the aircraft with its azimuth angle relative to runway 
centerline up to ^60°. The DME antenna which provides the range of the aircraft 
is generally co-located with the azimuth antenna. If the DME is located to the side 
of the runway, a simple transformation can still be used to obtain the aircraft’s 
position. The elevation antenna is located at the glidepath intercept point (GPIP), 
but is offset to the side of the runway; it provides the aircraft’s elevation angle 
up to 20° Thus, the aircraft has accurate position information in the volume of 
space within the limits mentioned above. 
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Consider a right handed coordinate frame with its origin at the phase center 
of the azimuth antenna, the x-axis along runway centerline and positive towards 
the runway, and the z-axis positive vertically upwards. Then if the position of 
the aircraft in this coordinate system (the MLS coordinate frame) is (Xq, yo> zq )> 
then the MLS signals would have the values given by the formulas below. 


R = 

• + 

2 , 2 1 
Yo + Zo 

(69) 

Az 

• —1 

-Yo 


= sin 

R 

(70) 



Zo 


El = 

= tan ^ 


( 71 ) 



V(jQ3-Xe)^ + (yo-yef ' 


where x^ and y^ are the x and y coordinates of the elevation antenna phase center 
in the MLS coordinate frame. The MLS signals received on the airplane are sim- 
ulated with additive white noise, bias and jitter. 

The onboard sensors used are three body-mounted accelerometers, attitude 
gyros for the pitch, roll and yaw angles, attitude rate gyros, barometric altitude 
and sink rate, airspeed and radar altitude. The accelerometers are mounted so 
as to measure the specific force along the aircraft's body axes. Let L^g be the 
matrix which transforms a vector from earth-fixed coordinates to body-axis 
coordinates, and L^g the matrix transformation from body-axis coordinates to 
the actual accelerometer axis; i. e. , it is assumed that the accelerometers are 
misaligned and are not exactly parallel to the aircraft's body-axes. Then the 
accelerometer reading, F, is modelled as 


F = Lab fB + ba + na , 

(72 ) 

^B - Lbb g , 

(73) 


where ag is the acceleration vector in body axes, fg the specific force vector in 
body axes, b^ the bias vector, n^ the noise vector and g is the gravity vector in 
the earth-fixed coordinate frame, given by 

= ( 0 0 32.2) . ( 74 ) 
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The earth-fixed coordinate frame is a right-handed coordinate frame with its 
origin at the GPIP on runway centerline, the x-axis positive toward the azimuth 
antenna, and the z-axis positive vertically downward. 

The attitude gyros for the pitch, roll and yaw angles are modelled as having 
an additive bias and additive white noise. The attitude rate gyros do not have a 
bias error but are modelled as being corrupted by additive white noise. The baro- 
metric altitude and radar altitude measurement errors are modelled by additive 
bias and additive white noise. On the other hand, the barometric sink rate and 
airspeed measurement errors are modelled by multiplicative noise. Thus, if 
Vjn is the measured airspeed, then 

Vm = ( 1 + e ) V , ( 75 ) 

where V is the aircraft’s airspeed and e is a zero mean white noise process. 

At high altitudes, the radar altitude accuracy is low due to non-linear effects; 
however, it is accurate at lower altitudes. On the other hand, note that since the 
elevation antenna is located at the GPIP, the nominal touchdown point of the air- 
craft is past the elevation antenna. Thus, during most of the flare maneuver and 
touchdown the elevation signal is not available to the aircraft. Hence, before the 
flare maneuver is started; the filter switches from using the elevation measure- 
ment to using the radar altimeter in order to obtain the position of the aircraft. 

Since at this point the altitude is low, the radar altitude provides reliable data. 

If a second elevation antenna located further down the runway is present, then this 
signal could also be used for flare and touchdown. 

The values for the standard deviation of the various noises and bias errors 
in the sensors were chosen to reflect current instrumentation standards [ 12 ] , 

C 13 ] ; the values are given in Table 1. 

B. FILTER EQUATIONS 

The measurements obtained from the on-board sensors and the MLS describe 
the motion of the aircraft and constitute the inputs to the filters. However, to use 
these measurements as the inputs to a Kalman filter, it is necessary to express 
them as linear combinations of the state variables. This can be done by pre- 
processing the measurements which are non-linear functions of the state variables. 
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Since the attitude of the aircraft, the attitude rates, the barometric and radar 
altitudes are linear functions of the state variables x and x{,, it is only necessary 
to pre-process the MLS signals, the accelerometer readings, the bare sink rate 
and airspeed. 

From the expressions for range, azimuth and elevation given in equations 
( 69 ) , ( 70 ), ( 71 ) , it is possible to solve for Xq , yo * z© . 


Yo = -R sin Az , 

(76 ) 

Xo = Xg sin^ El + a/x^ 0 sin^ El + ^ , 

( 77 ) 

Zo = tan El a/(Xo - Xef + (Yo - Ye)^ . 

(78 ) 

= (B? - Yo^) cds^ El - (^(Yo - Ye)^ + ) sin^ El . 

(79) 


Now, Xo , yo > Zq are linear combinations of the longitudinal and lateral state 
variables and x, and can be expressed as 

x'o = = -xts + + v^B ( 80 ) 

yo 

S'’” ' "ST ° '‘e * "s • ' ' 

Zo 

z'o = = -xts + , ( 82 ) 

where v^g , v^q , and va are additive white noise processes. Thus, these values of 
y’o z’o can be used as inputs to the filter and the expressions ( 80 ) , ( 81 ) 
and ( 82 ) can be used in the computation of the filter gains. 

The specific force measurements obtained from the body-mounted acceler- 
ometers can similarly be processed to obtain linear measurements. Using the 
relationships among the specific force, acceleration, and attitude rates it can be 
seen that 

Fx 

T? — ^BE ( 3 ) g ' %2 - u’ 

Uq 

= ej { Al XI + ui + nt w>t) + bax + Vax ( 83 ) 
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-^ + Lbe ( 2, 3 ) g - Tla - ^ + r 
T 

~ C 3 ( Ax + Bu + DOyy ^V ) + Xg + ^ 3 y ^3y 

F • 

+ Lbe (3, 3)g-11(j = -^-q 

= fj ( xt, + ut + nt, w^,) - xt4 + baz + Vaz 


( 84 ) 


(85) 


where g is 32. 2 /Uq , C{ is the i*^ column of the identity matrix, bgx , bgy and 
baz are the bias components in the processed versions of the specific forces in the 
X, y and z directions, respectively, and Vax . Vay and Vaz are vAiite noise com- 
ponents in the processed measurements. Thus, the left-hand-side of equations 
( 83 ) , ( 84 ) and ( 85 ) show the processing that needs to be done to the body- 
mounted accelerometer readings to obtain new pseudo measurements which can 
be expressed linearly in terms of the longitudinal and lateral state variables; the 
right-hand-side of the equations show the linear measurement model used. 

The barometric sink rate model was obtained directly from the expression 
given in equation ( 31 ) by normalizing the value by Uo and subtracting the predicted 
value of Tlz . The airspeed was linearized using 



(Of + a^,) 


3 


k 


1 


T (86 ) 

+ V , 


where v is a white noise process. 

The processed measurements pertaining to the longitudinal and lateral state 
variables can be separated and treated with the corresponding dynamical equations. 
Thus, the roll and yaw angles, the yaw rate and roll rate, the processed MLS y’o 
measurement and the processed body-mounted accelerometer Fy measurement 
were used in connection with the lateral equations of motion, the remaining 
measurements with the longitudinal equations. Treating each set separately, two 
filters were developed. The equations of motion were discretized as is shown in 
section IV, and the steady-state Kalman filter gains were computed. Then, the 
filtering equations were programmed in a non-linear aircraft dynamics simulation. 
The filter gains were further adjusted according to the closed loop performance 
obtained in the simulation runs. The form of the equations for the lateral filter is 
given below. 
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Vk = 

AAA 

^k " ^k ~ ^k » 

(87) 

A 

%+l = 

■ ^xk + Fuk ^k ♦ 

(88 ) 

A 

Wk+a 

II 

(89) 

% = 

A 

% ^x '’k ’ 

( 90) 

II 

A 

Wk + '^k • 

(91) 

^k+i “ 

A 

^k + ^ '^k > 

(92) 


AAA 

where , hy. are predicted values and 3^ , , Tj|^ are filtered values. A 

filter of the same form was used for the longitudinal dynamics. Various plots 
of the filter outputs obtained from the simulation runs are shown in section V. 
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IV. DEVELOPMENT OF THE CONTROL LAW 


The aircraft's motion was modelled by the longitudinal and lateral equations 
of motion as shown in section II. The model includes the response of the actuator 
servomechanisms and engine dynamics up to first order, and the effects of wind 
conditions on the motion of the aircraft. The wind model was developed in sec- 
tion n C using random gust, steady and shear wind components. The filtering 
equations were described in section m. In this section, these models are used 
to develop an automatic control law for the final approach till touchdown. The 
phases of flight considered are the localizer capture and track, steep glideslope 
capture and track, decrab and flare to touchdown. The lateral and longitudinal 
control laws are treated separately although some cross-coupling is present; i. e. 
some longitudinal variables are used for feedback in the control of lateral motion 
and vice versa. Digital design techniques are used in the development of the 
control laws by first discretizing the equations of motion into difference equations, 
and then using a quadratic cost criterion to obtain an optimal control law. The 
desired flight path is modelled so as to obtain satisfactory performance during 
the various phases of the final approach, 

A. DEVELOPMENT OF THE CONTROL EQUATIONS 

Consider the model developed for the lateral motion of the aircraft and the 
winds which affect this motion as given in equations ( 56 ) and ( 57 ). Note that 
the longitudinal equations of motion given in equations ( 37 ) and ( 68 ) are of the 
same form as the lateral ones. Thus, the development that follows applies to both 
models. The lateral equations are given below for convenience, 

X = Ax + Bu + DCw W + Tl , ( 56 ) 

W = Aw W + Bw “J • ( 57 ) 

These differential equations can be expressed equivalently in difference 
equation form provided that the control u remains constant over the sampling 
interval; for convenience, we shall further assume that Tlalso remains constant 
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over the sampling interval, although a more general case can be treated sim- 
ilarly. Hence, if 

u(t) = Uk , Tl(t) = % , tk^t<tk+i (93) 

where tk is the sampling instant equal to kT , T being the sampling period, 
then the differential equations in ( 56 ) , ( 57 ) can be integrated over the sam- 
pling period to obtain a discrete representation [ 14 ] , C 15 ] , [ ll] , 

^k+x = + r’Uk +Tw Wk + + u)ik • ( 94 ) 

Wk+i = % Wk + , ( 95 ) 

where xk and Wk are the values of x ( t ) and W ( t ) at the sampling instant tk, 
respectively, (u^k ^nd togk are zero mean white noise sequences with appropriate 
covariances, and 

$ = eAT = $( T ) , % = e^"'^ , ( 96 ) 


Jc 

> ds , r= , 

[T 

(97) 

rw = 

Jo ^( T-s) DC^ (s) ds . 

(98) 


Thus, the values of the state at the sampling instants can be obtained using 
the difference equations shown in { 94 ) and ( 95 ) . 

Depending on the phase of flight, the guidance law will require the aircraft 
to follow a certain flight path or flight condition. This desired flight path can be 
described by specifying the values of the state vector x corresponding to this 
flight condition. Thus, suppose that the vector z ( t ) describes the desired flight 
path; then the error in the actual flight path is 

e(t)=x(t)-z(t). (99) 

Since the objective is to follow the desired flight path, z, minimizing a cost function 
which is quadratic in the error, e, is appropriate; hence. 
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( 100 ) 


ptf 

J = I E Jo Ce' ( t ) Q e ( t ) + u' (t) R u ( t ) ] dt , 

where E is the statistical expectation operator. Now, the desired flight path can 
often be modelled by a differential equation. A model of the following form was 
used in the development here 


( 101 ) 
( 102 ) 


Z ^2 ^ ^2 * 

C ~ ?2 » 


where assumed to a gaussian white noise process. Thus, the desired flight 
path is modelled as a random process; i. e. , it is modelled as belonging to a family 
of trajectories which has given statistical properties. Inherent in this formulation 
is the fact that the desired flight path is not necessarily known a priori; i. e. , the 
future values of z are not necessarily known at the outset; hence, future values of 
the desired flight path are not needed to compute the current value of the control. 
Future values of the z would be necessary if it were modelled as a deterministic 
flight path [ 16 ] . Using this model for z, the equations of motion in ( 56 ) and 
( 57 ), the cost function, J, can be expressed in an equivalent discrete form [ 15 ] , 
[11 ]. 


Jj = I E J^e'k Qek + u'k RUk + 2 (e'k Ndk + e'k Muk + d’k SUk) 

Q = Jo $• ( s ) Q ( s ) ds , 

F , 

R = Jo r ( S ) Q r ( s ) ds + RT , 

M = JJ ( s ) Q r ( s ) ds , 

N “ Jo (®) Q^d (s) ds , 

F 

s = Jo r'd ( s ) Q r ( s ) ds . 


(103) 
( 104 ) 

(105) 

(106) 
( 107 ) 
( 108 ) 


The continuous and discrete cost functions J and Ji differ by a constant which does 
not depend on uk ; hence, the same control minimizes both cost functions. The 
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optimal control sequence which minimizes the cost function J given in ( 100 ) 
with the constraint that the control remain constant over the sampling period 
can be shown to be of the form [ 11 ] , [ 17 ] . 

) " Hw % - Hz Zk - ( Czk - % ) » ( ^09 ) 

where Xj^ , Wk and Tjj^ are filtered least mean square estimates. 

B. MODELLING OF THE DESIRED FLIGHT PATH 

The feedback control given in the last section requires that the model for the 
desired flight path, z, given in ( 101 ) and ( 102 ) be specified, since the values of 
z and Cz enter directly into the computation of the control commands. Note, 
however, that the matrix in the model is also used in the computation of the 
feedback gains, and thus must also be specified. 

The phases of the final approach considered here are the localizer capture, 
localizer track, glideslope capture and track, decrab and flare. The aircraft 
must perform these maneuvers in a constant airspeed, 3D mode. Thus, it is 
assumed that the aircraft is approaching the runway in level flight at a specified 
heading which intersects the runway centerline at a specified point. At some point 
on this flight path, the localizer capture mode is engaged, and the automatic con- 
trol law controls the aircraft through the various phases of the final approach till 
touchdown. 

The localizer capture engage logic was developed so that the initial rolling 
tendency of the aircraft is to bank away from runway centerline, thus avoiding 
overshoot tendencies. Now, let ^ be the yaw angle estimate and y* the estimate 
of the normalized distance of the aircraft from runway centerline. Now, assuming 
that the remaining variables in equation ( 109 ) are small, the aileron command, 

6A, will be given by 

6a = -He ( 1, 2 ) ij; - He ( 1, 6 ) y» . ( 110 ) 

Thus, localizer capture should be engaged when the aircraft crosses the curve ob- 
tained by setting 6 a equal to zero. Further observation of the situation leads to 
the following criterion. 
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( 111 ) 




He (1,2) 

y 

^Uo 

'ij 

He(l,6) 


where Uo is the inertial speed of the aircraft, and y is the distance from runway 
centerline. Thus, localizer capture is engaged when the criterion in ( 111 ) is 
met. It should be noted that this criterion is useful as long as y has a small error 
as would be the case when the MLS signal is valid; if the position estimate has large 
errors, capture could be initiated too early or too late to reduce overshoots. 

The desired flight path for the lateral motion of the aircraft was chosen to 
be the runway centerline, with the desired cross-runway velocity being zero. 

Under no wind conditions, the desired yaw and roll angles as well as their rates 
of change are zero; however, if a cross wind is present the aircraft should crab 
into the wind and then decrab before landing while keeping a small aerod 5 mamic 
sideslip angle throughout the flight. It is also desired that any overshoot or under- 
shoot during localizer capture be as small as consistent with the other requirements. 
The desired flight path in equations ( 101 ), ( 102 ) was modelled using 

Az = A , Cz = z - Az , ( 112 ) 

for the lateral dynamics. The discrete form was used to update the desired flight 
path. 


Kk+i = ^ zk ^zk • ( ) 

In order to have the aircraft crab into the wind rather than use a bank angle, 
the yaw angle weight was set to zero and the desired inertial sideslip angle was set 
equal to the negative value of the estimate of the sideslip angle due to steady cross- 
winds. 


Z]j3 - - ^ 

To reduce overshoots, the closed loop damping on lateral offset was obtained by 
using a non-zero weight on the cross-runway velocity. Decrab was obtained by 
setting the desired roll rate to be proportional to the deviation in the yaw angle 
from the runway heading, and, simultaneously, setting the desired yaw rate pro- 
portional to the yaw angle deviation, when the aircraft crosses the decrab altitude. 
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An exponential ”easy-on” was used to introduce the decrab commands in order 
to avoid sudden changes in the commands and the ensuing transient response. 

The desired flight path for the longitudinal variables, was specified 
in a similar manner. To obtain a constant airspeed of Uo, the normalized de- 
sired inertial speed along the stability x-axis was set equal to the difference 
between the normalized desired airspeed and the normalized wind velocity esti- 
mate along the stability x-axis. 

z^ = 1 - €’ Wi , ( 115 ) 

To follow a glideslope of angle y© » the inertial velocity must point along the glide- 
slope; the desired angle of attack of the inertial velocity was set so as to obtain 
this proportion between the x^ and Zq components of the inertial velocity* 

+ tan Yq L^i ) ( 1 + Z 3 ) +(^£32 + tan Yq Lisy / 

( f <33 tan Yo- , ( 117 ) 

where L^j is the estimate of ( i, j ) and 3 is the estimate of the inertial sideslip 
angle. 

To obtain a 3D control law the longitudinal model developed in section II was 
slightly modified; however, the filter development is based on the unmodified model. 

To make the control law independent of the x position of the aircraft, this variable 
was excluded from the control model; thus, the number of state variables were 
reduced by one. On the other hand, the altitude variable (xte) was modified to 
represent the altitude error from the glideslope. This error was defined as the 
difference between the aircraft’s actual altitude and the altitude of the point on the 
glideslope corresponding to the aircraft’s actual x position; i. e. , 

He = Z’e + tanyo X’e , ( 118 ) 

where X’e and Z’^ are the x and z components of the aircraft’s position in the earth- 
fixed coordinate frame normalized by Uq. After some manipulation, a differential 
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equation for x^e in terms of the other state variables can be obtained. 
xie = -xti + secyo xta + % 

T\le = ( I<3i + tan yo Lu ) ( 1 + xts ) + ( Lq3 + tan yo Li 3 ) xt3 + ( L 33 + tan yo Li 2 ) 3 

+ X(,i - secyo X(j , ( 119 ) 

where Ljj is L^s ( i, j ) , and Pis the inertial sideslip angle, , from the lateral 
model. The model A 2 ^,used was the same as the one given in [ 11] with the change 
that was omitted and z^s was given by 

= -z^i + sec yo zi 3 + • ( 120 ) 

To obtain Cz-te . z^s was first computed by 

zie = Z'ed + tanyo , ( 121 ) 

where Z'g^ is the desired normalized sink rate, and is the estimate of the air- 
craft’s inertial velocity along Xe. Thus, z^e is obtained; then is computed by 
solving equation ( 120 ) . 

The control feedback gains for the longitudinal control law were computed 

changes in the 
simulation runs. 


using the modified model for the motion of the aircraft. Various 
gain values were made to improve performance as obtained from 
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V. RESULTS 


The performance of the automatic landing system described in the preceding 
sections was evaluated using a non-linear digital simulation of a Boeing 737 air- 
craft. Since the automatic landing control law was developed using linear methods 
of modem control theory and was tested with a linear simulation of the aircraft 
using perturbation equations, evaluating its performance with a more accurate, 
six-degree-of- freedom, non-linear simulation of the aircraft was considered to 
be a realistic approach to its performance evaluation. The sensor models used 
in this simulation were also more accurate in their generation of error signals 
and noisy readings. 

The flight paths considered in the evaluation consist of steep final approach 
paths starting at localizer capture and ending at touchdown. Thus, at the beginning 
of the simulation, the aircraft is assumed to be in level flight, approaching the 
runway at approximately 61. 73 m/sec (120 knots) on a straight line which intercepts 
the runway centerline at a specified track angle and a corresponding yaw angle. 

When the criterion given by (111) is met, the localizer capture mode is engaged and 
the control law starts to align the aircraft with the runway. Similarly, when the 
glideslope capture criterion is met, the aircraft starts the glideslope capture man- 
euver, which brings the aircraft to a steep glideslope (6®) by pitching down and 
increasing the sink rate to its appropriate value depending on wind conditions. 

Since the control law tries to maintain a constant airspeed, in a head wind condition, 
the ground speed will be lower, and to stay on the glideslope, the aircraft will have 
a lower sink rate; in a tail wind condition, the situation would be reversed resulting 
in a higher sink rate. 

By appropriate choice of the initial altitude, the localizer and glideslope capture 
modes can be initiated at the same time. In this case, both capture maneuvers 
occur simultaneously; and the aircraft follows a curved path in the horizontal plane 
while also following a curved path in the vertical plane. Thus, the aircraft banks 
away from the runway to reduce its cross-runway velocity while it is pitching down 
to increase its sink rate. This simultaneous capture of localizer and glideslope is 
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desirable as it reduces the amount of time (and space or distance) used in the 
capture maneuvers. It allows aircraft to perform the capture maneuver closer 
to the runway, thus using less space on runway centerline which now becomes 
available for other aircraft. Thus, the aircraft can perform close-in simul- 
taneous localizer and glideslope captures. The simulation runs made for the 
evaluation of the control law were chosen so as to perform simultaneous captures. 

Another factor which affects close-in captures is the amount of time it takes 
for the aircraft to arrive at a stable condition on the glideslope and localizer, 
i.e. the capture settling time. Thus, if the settling time is small the aircraft 
is ready to start the subsequent maneuvers, in this case decrab and flare. 

This would allow the aircraft to start the capture maneuvers closer to the 
runway. Thus, in the development of the control law, an effort was made to 
complete the localizer and glideslope capture in a short period of time by avoid- 
ing large overshoots and undershoots, and large time constants for the closed- 
loop system. As can be seen from the plots of the various parameters obtained 
from the simulation runs, it takes approximately 40 seconds for the aircraft to 
complete the simultaneous localizer and glideslope captures and settle in a steady 
flight condition on the glideslope, starting from level flight with localizer intercept 
angles of up to 40°. 

In cases where the wind velocity has a cross-wind component, the lateral 
control law automatically crabs the aircraft into the wind, while remaining on 
runway centerline; thus, the aerodynamic sideslip angle is kept small. This is 
achieved by estimating the magnitude of the cross-wind component and feeding this 
information back through the control commands. The various plots obtained from 
simulation runs where a cross-wind is present show that the control law achieves 
the correct crab angle for the various wind conditions considered as can be seen 
by checking the yaw angle of the aircraft relative to runway centerline and its 
■offset from runway centerline. Also, note that no excessive overshoots are pre- 
sent in localizer capture in these cases, irrespective of the direction of the cross- 
wind. This performance is also due to the estimation and feedback of the wind 
velocity. 
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When the aircraft has a crab angle due to cross-winds, the decrab mode is 
engaged when the estimated altitude reaches a specified altitude. This decrab 
altitude is currently set at 76.2 m since the aircraft is on a steep glideslope and 
its sink rate is higher than the eonventional 2. 5° - 3° glideslopes. The decrab 
law is designed to align the aircraft's x body axis with the runway centerline; 
i. e. it reduces the crab angle by rolling into the wind while the yaw angle is 
reduced to zero. This is achieved by commanding the aircraft to roll into the 
appropriate side until the yaw angle reaches zero. A yaw rate command pro- 
portional to the yaw angle is also used to aid in the maneuver. The simulation 
runs conducted show that the aircraft reduces its yaw angle and rolls into the 
wind while remaining on runway centerline when the decrab mode is engaged; 
the results under cross-winds accompanied by gusts also appears satisfactory 
in the simulation runs made. 

As in the case of the decrab altitude, flare is initiated at a higher altitude 
than it would be for shallower glideslope angles due to the fact that the sink rate 
is higher with steep glideslopes. Thus, the flare mode is engaged when the air- 
craft's estimated altitude reaches 50. 50 m . Since the sink rate is higher on a 
steep glideslope and flare is a flight critical phase final approach, initiating flare 
at a higher altitude has the advantage of lowering the sink rate to more conven- 
tional levels at a higher altitude. This provides the advantages of a slower and 
more stable flare maneuver, and leaves more time and altitude for the pilot to 
monitor the progress of the maneuver with ease. In this mode, the aircraft 
follows a flare path with a specified vertical profile until touchdown. The flare 
path is generated on line with the independent variable being the distance to the 
GPIP. The simulation runs show that the flare law achieves satisfactory landings 
under the various wind conditions considered. 

The digital nonlinear simulation ALERT was used for the aircraft dynamics 
of the B-737 and the various sensor dynamics and errors. The aircraft dynamics 
model in ALERT uses the six-degree-of-freedom nonlinear equations of motion. 
The aerodynamic forces are generated using the stability derivatives of the B-737 
aircraft for an airspeed of 61. 73 m/sec (120 Imots). A detailed description of the 
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simulation can bo found in [ 18 ] . The basic update rate for the simulation of the 
aicrcraft dynamics was 20 samples per second, corresponding to a period of 50 
msec. . The digital control law operates at a sampling rate of 10 Hz, or twice as 
slow as the basic simulation update rate. Thus, the sensor measurements were 
sampled every 100 msec, although they were available every 50 msec, in the 
simulation. 

The MLS signals of azimuth, elevation and range, were simulated by cor- 
rupting the true values by an error signal. The error signal contained noise 
and bias to model the effects of calibration errors and misalignments in the 
installment of the antennae. The jitter that has been added to the MLS signals 
to avoid interference was also modelled and simulated in the error signal. The 
body-mounted accelerometers were modelled so as to include the effects of 
misalignment relative to the aircraft's body axes, errors in scaling, bias and 
noise. Thus, the true values of the specific force components along the body 
axes are rotated to obtain the readings of the misaligned specific force compon- 
ents, then scaling error effects are introduced before the bias and noise are 
added. The attitude rate gyros and airspeed indicator were modelled as being 
corrupted by noise, whereas the attitude gyro, barometric altitude and sink rate 
and radar altimeter measurements contained bias errors and noise; the noises 
in airspeed and barometric sink rate were introduced as multiplicative noise. 

The values for the standard deviation of the noises and other errors reflect cur- 
rent instrumentation standards [ 12 ] , [ 13 ] . Table 1 shows the values of the 
various parameters used in the error models of the sensors. 

Numerous simulation runs were conducted using the ALERT simulation of the 
aircraft and navigation aids in partial evaluation of the overall performance of the 
control law developed. Various parameters describing the response of the closed 
loop system (aircraft dynamics, sensors and control law) were plotted. Figures 
5-20 show the plots of these parameters for some of the simulation runs con- 
ducted. These simulations show the response of the aircraft under steady winds 
of 5. 14 m/sec (10 knot) velocity from various directions: head winds, tail winds. 
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cross winds and other combinations, turbulence of various intensities with and with- 
out steady winds, and under wind shear conditions. For the purposes of this inves- 
tigation, wind gusts with a standard deviation of .61 m/sec in the forward and 
lateral directions ( x and y body axes) and .30 m/sec in the vertical direction were 
denoted as average turbulence levels; doubling the standard deviation of the wind 
gust velocities in all three directions was accepted as a high turbulence level. The 
wind shear conditions were simulated starting at an altitude of 60. 96 m and ending 
at touchdown. From 60.96 m to 30.48 m of altitude the steady wind is decreased 
linearly at a rate of 2.06 m/sec (4 knots) per 30.48 m (100 feet); from 30.48 m to 
touchdown the steady wind changes at a rate of 4. 12 m/sec (8 knots) per 30. 48 m. 
Thus, the steady wind component changes 6. 17 m/sec (12 knots) in 60. 96 m pro- 
ceeding touchdown with the rate of change increasing as the aircraft gets closer to 
touchdown and during flare. 

Various of the simulation runs made under combinations of steady winds, 
turbulence levels, turbulence patterns or histories and wind shear conditions are 
shown in Figures 5 - 20 . Localizer captures with differing intercept angles 
were also simulated to examine the capture response under differing initial con- 
ditions. Simulations of the final approach with nominal airspeeds of 61. 73 and 
64. 30 m/sec were made, with the majority of the runs being at the 61. 73 m/sec 
design airspeed. On various occasions, the initial conditions were offset from the 
desired, nominal or trim values; the initial values of the filter estimates were also 
intentionally set to incorrect values to assess their effects on performance. Various 
other parameters were also changed during the evaluation to check their sensitivity 
and to improve the response. 

From the simulation runs made, it is seen that the overall performance of the 
control law during the final approach, starting from localizer and glideslope capture 
to touchdown, is satisfactory under the various wind conditions considered. The 
simultaneous localizer and glideslope capture maneuvers are performed and com- 
pleted in a short period of time and without excessive overshoots in the presence 
of head winds, cross winds and turbulence as well as under no wind cases. When 
the aircraft is approaching the runway centerline at a 30^ track angle localizer 
capture is initiated at 361 m offset from centerline; obseirvation of the plots for 
this initial condition (as well as others) shows that localizer capture is completed 
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in less than 40 seconds, so that for all practical purposes the aircraft has set- 
tled on runway centerline. It can further be seen that the lateral overshoot in 
these cases is small, sometimes barely noticeable, and rarely exceeding 7. 6 m 
at its peak. The lack of an appreciable overshoot in lateral offset can be traced 
to the use of wind velocity estimates in the feedback control; when the initial 
cross wind estimate (at capture initiation) has a 5. 14 m/sec error, larger over- 
shoots approaching 15 m are obtained before the error in the estimate is reduced 
to lower levels. These situations can be avoided by running the filter open loop 
before the control law is engaged. On the other hand, a more basic approach 
would be to use a sensor measuring the aerodynamic sideslip ai^le (such as a 
P - vane or other sensor) to reduce the time constant related to the cross-wind 
estimate. This would reduce the error in the wind estimate in a shorter period 
of time and reduce any overshoot as well as aid in gust alleviation and wind shear 
accomodation. 

Similar remarks can be made about the glideslope capture which is performed 
at the same time as localizer capture. From the various plots, it can be seen that 
the aircraft settles on the glideslope in about 20 seconds as the aircraft is more 
responsive vertically and the offsets involved are smaller. It may be noted that, 
imder no wind conditions, the aircraft has a vertical offset which gets smaller as 
it gets closer to the runway. This is due to the bias or misalignment simulated in 
the MLS elevation signal which is not accounted for by the filter. A similar offset 
with a decreasing trend can be seen in the lateral direction due to the bias in the 
MLS azimuth signal simulation. The effects of these offsets are small, predictable, 

and are considered of minor significance relative to other offsets. 

The overall performance during flare and touchdown under various wind condi- 
tions is shown in Figures 5 - 20. The sink rate which holds a constant level on the 
glideslope gradually decreases to a smaller value and tries to maintain this lower 
rate until touchdown. The desired sink rate at touchdown was chosen to be . 64 m/ sec 
(2. 1 ft/sec) at 61. 73 m/sec. This constant sink rate command until touchdown was 
used to reduce the touchdown dispersion when wind conditions tend to make the air- 
craft float past its nominal touchdown point. The plots of the aircraft’s sink rate 
under various wind conditions show a satisfactory response during the various phases 
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of flight and flare. In the simulation runs made, the touchdown sink rate has a mean 
or average value of . 6C m/sec and a standard deviation of . 26 m/sec. The average 
value of the touchdown points was 416 m from the GPIP with a standard deviation of 
34. 5 m. Table 2 shows the touchdown sink rate and the deviation in the touchdown 
point for some of the simulation runs made under differing wind conditions. It 
should also be noted that the other parameters relevant to the touchdown condition 
such as pitch attitude, yaw angle, offset from runway centerline, etc. also had 
satisfactory values as can be seen from Figures 5 - 20. 

In comparison to simulations of the autoland control law currently used with 
the third order complementary filter or the time varying Kalman filter C 18 ] , 

C 19] , [ 20 ] , the digital control law described above has smaller overshoot 
peaks and requires less time to settle to a steady condition during the localizer 
and glideslope capture maneuvers. It should also be noted that the latter law per- 
forms the localizer and glideslope captures simultaneously and follows a steep 
glideslope. The latter law computes estimates of the wind velocity components 
and uses these estimates in the feedback control commands whereas the current 
autoland law uses a wind shear detector. A difference in the design approach is 
due to the fact that the current autoland control law was developed as a continuous 
or analog system and was implemented on a digital computer using a 20 Hz sam- 
pling rate, whereas the latter control law was developed as a digital or discrete 
system based on a 10 Hz sampling rate using modem digital control techniques . 

In conclusion, the digital automatic landing guidance and control law, as seen 
from the simulation evaluation described above, has satisfactory performance in 
the final approach to landing phases of flight. 
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VI. SUMMARY 


A digital guidance and control law for a steep final approach was developed 
using modem digital control techniques for operation in a 3D constant airspeed 
mode. To take advantage of the high precision MLS guidance information, a con- 
stant gain digital Kalman filter which combines the MLS information with various 
on-board sensor outputs was developed to obtain accurate estimates of the air- 
craft’s position, velocity, attitude and the wind velocity components. The wind 
velocity estimates are used in the feedback control law along with the other dy- 
namical variables for gust alleviation and wind shear detection and control. The 
control law considers simultaneous localizer and glideslope capture, localizer 
and glideslope tracking, crab and decrab, and flare to touchdown. The digital 
filter and control law were simulated on a digital computer with a non-linear 
simulation of the B-737 aircraft, on-board sensors and the Microwave Landing 
System. Simulation runs made under varying wind conditions show satisfactory 
performance of the automatic control system throughout the final approach. 
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Figure 6 Final Approach Simulation: 20^ localizer intercept 
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Figure 7 Final Approach Simulation: 30° localizer intercept 






















Figure 11 Final Approach Simulation: tail wind, average turbulence 
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Figure 20 Final Approach Simulation: 64. 3 m/sec (125 knot) airspeed, average gusts 


TABLE 1 Sensor error model parameters 


Sensor 

Noise standard 
deviation 

Bias 

Comments 

Attitude gyro 

.229° 

.229° 


Rate gyro 

. 02®/sec 

0 


Body-mounted 

accelerometers 

1% of g 

1% of g 

. 229° misalignment 
and 

. 25% scaling errors 
also modeled 

MLS - azimuth 

o 

O 

. 0125° 

jitter added 

MLS “ elevation 

o 

o 

. 0125° 

jitter added 

MLS - DME 

2.29 m 

2.29 m 


Radar altimeter 

.305 m 

. 305 m 


Barometric altimeter 

. 305 m 

1. 52 m 


Barometric sink rate 

. 305 m/sec 

.61 m/sec 


Airspeed indicator 

2% 

0 

noise is multiplicative 


g = 9. 81 m/sec^ 


58 








TABLE 2 Touchdown Statistics 


Wind Conditions 

Touchdo\vn Conditions 

head wind (+) 
tail wind (-) 

cross-wind 

turbulence 

shear 

sink rate 
(m/ sec) 

deviation from 
average touch- 
down pt. (m) 

0 

0 

0 

F 

.75 

21.7 

0 

0 

average 

F 

.91 

25.3 

0 

0 

high 

F 

1.08 

31.4 

4 - 

0 

0 

F 

.29 

-4.5 

+ 

0 

average 

F 

.37 

-32.6 

+ 

0 

high 

F 

.40 

-64.6 


0 

0 

F 

.87 

-21.0 

- 

0 

average 

F 

.43 

-9.1 

- 

0 

high 

F 

.56 

67.7 

0 

- 

0 

F 

.84 

-8.5 

0 

- 

average 

F 

.60 

2.19 

0 

” 1 

high 

F 

.65 

38.5 

0 

1 

+ 1 

0 

F 

.39 

-18.2 

0 

+ 

average 

F 

.75 

47.0 

0 

+ 

high 

F 

1.27 

54.9 

+ 

- 

0 

F 

.83 

1 

00 

CO 

+ 

- 

average 

F 

.73 

-26.2 

+ 

- 

high 

F 

i 

.80 

-38.7 

+ 

0 

0 

T 

.31 

-31.7 

- 

0 

0 

T 

.46 

J 

-12.4 


average sink rate: -665 m/sec (2.18 ft/sec) 
sink rate standard deviation: .26 m/sec (.85 ft/sec) 
average touchdown point: 415. 7 m (1363. 8 ft) 
standard deviation of touchdown point: 34. 51 m (113. 2 ft) 

All steady winds have a magnitude of 5.14 m/sec (10 knots) 
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APPENDIX 


The aircraft equations of motion given in section n are in state variable 
form. The relation of the elements of the matrices used in these equations to 
the stability derivatives of the aircraft are given here. The lateral equations 


repeated here for convenience. 


X = Ax + 

Bu - 

Dw + T] , 

( 56 ) 

< 

II 

W 

+ tu , 

• 

( 57 ) 

w = Cw 

W 

+ § . 

(55) 


The elements of the matrices in these equations are given below in terms of the 
aircraft’s stability derivatives, and other aircraft parameters. 

" q^Sb ’ 

= fi fg - fs , fe = 

where S is the wing area, b the wing span, qa the dynamic pressure, Ijcx» Ixz 
Igz the aircraft’s moments of inertia in stability axes, m the aircraft's mass 


Aj4 - 1 , 

Ai5 ~ 

tanBo 

9 Aas - 

secBo , 

. g COS0Q 

‘ u„ 

9 

A33 = 

^yP 

h 

A - ^yp 
A34 - ^ 

Cyr 

Aas - ^ 1 

9 

A37 “ 

Cy6 R 


A43 = 

+ fs 

( Cn6 

+ CnTP)] . 


A44 ” ( C^p 

fs + f 3 

^np ) 

9 



^zz j _ ^xz 

Qo Sb ’ ® Qo Sb 

p\Uq 

QoS 
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A 46 = ^ ^3 Cnr ) * A47 = ( Ci^Yl ^2 + h Cn6R) * 

Ab 3 = ( Cn 3 ■*■ CnT3 ^ ^^3 ) ♦ A54 = ( C^p + fa A44 ) , 

A55 “ ( Cyiy fa A45 ) , A57 = ( Cn 6 R + ^ A47 ) , 

Agi = -sino^j , Ag 2 = cosYo » Aas = 1 • 

The elements of the matrix A for which an expression is not given above are 
zero; similarly, only non-zero elements of the remaining matrices will be given 
below. 

Bai = ( Cy 6 A ^spa ^y^sp ) » 

^41 “ |^( fg + £3 Cn6A ^ Cspa ( ^^-sp fa + fa Cnfep 1 

B51 = [Cn6A ( f^-t6A fs fs 

^ spa ^f^n6sp ^ ( C-t,6sp fa Cn6sp)^J * 

B73 = 1 . 

Ebi = A33 , Eba = A34 > D33 = Ass » 

D41 *= A43 , D42 = A44 t D43 = A46 9 

Ebi “ A53 , Efes = A54 9 Efea = A55 

The parameters related to the lateral wind model are given below. 

Awi 3 ” 1 » 

A ^ A - - 2Uo_ 

^W2i - -( t Aw 33 - 
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A - _ "Uo 


II 

1 

Q 

< 

3 1 

n Uo Lv 

Bw 3 i - av ( 2 v/P - 1 J ' . 

* HL. V 

M 

II 

/jJUo ' 

^ 3 Ly 

R. - '^WTT 

Bw 42 

y 8 ^ ttLw^i ' 

^ Lw ^ 4 b 

B\v 53 == 1 » 

®W 54 == 1 

^wll " 1 > 

C ^15 =1 , 

C^y 23 “ sinc^Q 

» ^W 24 ~ COSCYO 

Cw 33 = coscyo 

> Cw 34 = -sincvo 


The aircraft's longitudinal equations of motion were given in ( 37 ), and are 
repeated here along with the wind model for the longitudinal components of the 
wind velocities 

(37) 
(68 ) 

The elements of the matrices involved are given below in terms of the 
stability derivatives and the following intermediate parameters. 


X(, = Atxi + ut + Dt + lU, 

W4 = Awl Wl + U , 'Ni = Cwl Wt • 


Aw32 

Aw44 

Awsa 


TTUo 

3b 

TTUo 

4b 


= 1 
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cva = 


mg BinBp 
qoS 


»4 = 1/ (Ota + Ci^ ) , 

qo = I P Uo^ . 

where S is the wing area, c the mean aerodynamics chord and p is the air density. 

Using these variables, the matrix elements are given below. 

aisi 

at3i =0/2 04 , , ai2S = -«4(Clu + 2Clo ) . = - “4 (CLo- + Cr>o 1 . 


= 1 


a, 1 

- — . at® “ ( -Cdu - 2 Cdo ^Txu 2 Ctxo ) 


Dor Lo 


ats7 - 


Tx 

O3 


afee = 


^D6s 

03 


Ol = 


mg cosoq 
qoS 


o’a = 


mUp 
qo S 


Of s “ ^yy/ qo 


Uoi = «4 ( “3 - CLq ) , at37 = «4 Ctz . atas = - 0f4 Cl6s 

at4i = ®4 Cmof 02 /os , at48 ~ 7^ ( Cj^u + 2 Cmo + Cj^q; atas ) » 


at43 = 
at47 = 


( ^ma ^mTor ^ma ) ♦ a^44 - ( C^iq + CjjjQr ^-t34 ) • 

^mT + Cmor ®-t37 _ Cj^gs + Cmor 

^ 


^t 77 “ • 5 


3<^7 9 — *2985 




* C f, + C • 

Doe , ^ , moe ma <31 

t b^x - - a4 CL6e * b^x ~ — 


<^3 






dtgg “ a^as ♦ 


“ 3-t33 


dt32 - at33 


d^33 - 0^4 ( - Cl/J ) 9 
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“ ^-^4 2 * dU3 “ ^43 


- 


^mq + ( ^33 ” ^ ^ ^mb 


P^B 


A'W't/lS 

Aw-tsi 

Aw'tsx 

'^W'C44 


-(^f 

IjW 


ttUq 

4b 


ik 

Lu 


Aw-t32 


2 Up 

Lyv 


■^w-taa 


-TTUq 

4b 


Aw-k? “ 1 


V3 Up 
Lw 


Bwisi = (1-V12 ) (-^)^ , 

Lw 

Bw-tTS = 1-° * 

Ln/, 

Bw-te* ~ ^ ’ B^^73 = 1 


C^tii = -sino-o , 

Cwtis ~ -cosYo » 

Cwt3l = -COSQ'o , 

Cw-tss “ ~sinY) i 

C . _ _ ^Uo 


— COSQJq 


Cw-Li 6 “ sinYo 


Cwt34 = sino-o 


Cwtss = -cosYo 


r , - "Bo 

^W^33 - 
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